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Abstract - Nonlinear wave patterns generated by the flow of polariton condensate past an ob- 
stacle are studied for quasi-one-dimensional microcavity geometry. It is shown that pumping and 
nonlinear damping play a crucial role in this process leading to sharp differences in subsonic and 
supersonic regimes. Subsonic flows result in a smooth disturbance of the equilibrium condensate 
around the obstacle whereas supersonic flow generates a dispersive shock wave in the flow up- 
stream the obstacle and a long smooth downstream tail. Main characteristics of the wave pattern 
are calculated analytically and analytical results are in excellent agreement with the results of 
numerical simulations. The conditions for existence of stationary wave patterns are determined 
numerically. 



Introduction. — Recent experimental progresses in 
studying the microcavity polaritons have lead to a huge 
growth of interest in their collective dynamics (see, e.g., re- 
view articles [l][2] and references therein). Polaritons pos- 
sess an extremely small effective mass which allows their 
condensation at temperatures much greater than that of 
ultracold atomic vapors. Besides that, parameters of the 
polariton superfluid can be easily tuned with the use of 
resonant lasers. However, polaritons have a finite lifetime, 
and to maintain their steady-state population a continu- 
ous pumping is required. Experimentally, it is observed 
3,4 that above a threshold pumping strength an accu- 
mulation of low energy polaritons is accompanied by a 
significant increase of spatial coherence that extends over 
the entire cloud of polaritons which can then be described 
by a single order parameter (polariton condensate wave 
function) obeying an effective Gross-Pitaevskii equation. 
On the contrary to the atomic condensate situation, the 
density of the polariton condensate is not an arbitrary 
parameter anymore. Instead, it is determined by the con- 
dition of balance between pumping and dissipation pro- 



cesses. This implies that the non-conservative effects can 
play a crucial role in the condensate's nonlinear dynam- 
ics. For example, the generation of oblique solitons by the 
flow of a condensate past a localized obstacle has been 
observed [5] at subsonic speed which is impossible in the 
case of conservative atomic condensate. The formation of 
oblique solitons followed by their decay into vortex streets 
in a non-uniform cloud of polariton condensate has been 
studied in jfjj. Other geometries are also of great inter- 
est. In particular, quasi-lD flow of polariton condensate 
along "quantum wires" was studied experimentally in [7]. 
In the atomic condensate case such a flow past an obstacle 
leads to generation of dispersive shock waves propagating 
upstream and downstream from the obstacle [8 10 , as it 
was observed in the experiment 11 and explained theo- 
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retically in 12 . However, it is easy to see that this theory 
cannot be applied to the polariton condensate flow which 
density must be fixed (if pumping and dissipation are bal- 
anced) far enough from the obstacle thus preventing for- 
mation of jumps in the flow parameters. This means that 
in the dissipative case the dispersive shock waves gener- 
ated by the flow must always be attached to the obstacle 
and relax to the steady-state flow far enough from the ob- 
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stacle. This qualitative difference between wave patterns 
in conservative and dissipative cases makes the theory of 
dissipative flow much more complicated and this Letter is 
devoted to its development. 

Theoretical model. — Several models have been sug- 
gested for theoretical description of polariton condensate 
(see, 



e.g. 



13 



They were based on various generaliza- 
tions of the Gross-Pitaevskii (GP) equation. Here we have 
to take into account the effects of losses and pumping of 
polaritons on generation of dispersive shock waves by the 
flow of condensate past an obstacle. To this end, we will 



use the simple model introduced in 14 where nonresonant 



pumping due to stimulated scattering of polaritons into 
the condensate and their linear losses are described by the 
effective "gain" term d t i/> = 7^ where ip is the polariton 
condensate wave function. The saturation of gain is mod- 
eled by the nonlinear term dftp = —T\ip\ ip which brings 
the condensate density into equilibrium with the one of 
external reservoir, p = \ip\ 2 — 7/r. The dispersion of the 
lower branch of polaritons in the effective mass approxi- 
mation and the nonlinear interaction due to the exciton 
component of polaritons lead to the following generalized 
GP equation 

ifl* + \i>xx - l^l 2 ^ = v(x)i> + *( 7 - r|^| 2 )^, (1) 

(written in standard non-dimensional units) where V(x) 
denotes the potential of the obstacle and x is the direction 
of the flow of the condensate. The theory developed below 
can be generalized to other forms of nonlinear gain as, 
for example, the model considered in 15 , and the results 
remain qualitatively similar. However, to be definite, we 
shall consider here the model for which the description 
looks especially simple. 

It is easy to see that eq. ([I]) without potential (V(x) = 0) 
admits a plane wave solution which constant amplitude is 
determined by the gain and damping coefficients, 



fi = 7/r + \u, 



(2) 



where uq denotes the uniform flow velocity of the conden- 
sate and po = 7/r is its density. The analysis of modula- 
tion stability of such plane waves with respect to harmonic 
perturbations cx exp[i(Kx — fit)] shows that the distur- 
bance propagates along the wave ^ with the dispersion 
law 

fl = u K - i-y ± y/K 2 {a 2 + K 2 /4) - i 1 . (3) 

The expression under the square root vanishes at the wave 
number 

K c = a [2(^1 + T 2 - 1)1 (4) 

which separates two different regimes of evolution of a 
harmonic perturbation; see the plots in fig. [T] In the 
conservative limit 7, T — > 0,7/r — const we reproduce 
the standard Bogoliubov dispersion law with sound veloc- 
ity c s = a = y/pTJ- We shall call this parameter "sound 
velocity" also for small 7, T <IC 1; it corresponds to the 




Fig. 1: (Color online) The real and imaginary parts of Q as 
functions of K (see eq. (|3|) are shown in the case 7 = T = 0.05 
and tto = 0. 



almost linear part £1 ~ c s K of the dispersion law for 
K c k, Ta <C K <C a (see fig. [I]). Notice that, as long 
as repulsive nonlinearity is considered, the plane wave so- 
lution ([2]) is modulationally stable, since Im(O) < for 
any ii'-value. 

It is convenient to transform eq. |l]) into a hydrody- 
namic form by means of the substitution 



ip = y/p exp I i / u(x',t)dx' 



-ifit 



which yields 



Pt + (pu) x = 2Tp(p Q - p), 



. Px Pxx 



= -V x 



(5) 

(6) 
(7) 



The last term in the left-hand side of ((7| describes the 
dispersion effects; the right-hand side of Ml describes the 
gain and loss effects in the system; the other terms in 
eqs. ([6]),([7| have standard hydrodynamic meaning. 

Hydraulic approximation. In what follows we 
assume that the size I of the obstacle is much greater 
than the "healing length" (equal to unity in our non- 
dimensional variables). In our numerical simulations we 
will use an obstacle potential in the form 



V(x) = V m ex V (~x 2 /l 2 ), Z>1. 



(8) 



All numerical plots below are obtained for 1 — 5. We 
are interested in stationary patterns with p t = 0, Ut = 
generated by the steady flow of the condensate past a 
localized obstacle (V(x) — 5- as |x| — 5- 00), i.e. p and u 
must satisfy the boundary conditions 

p —> po, u — > u as I as j — > 00. (9) 

Then eq. ^ can be reduced to 

(pu) x = 2Tp( Po - p) (10) 

and the stationary eq. ^ can be integrated once to give 



1 2 1 I P X 

2 U + P + ^. 



Ap 



+ V{x) = \u 2 + Po (11) 
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For I 3> 1 it is natural to assume that the wave pattern has 
the characteristic length about I and, hence, the dispersive 



terms in (111 having higher order derivatives are negligibly 



small compared with other terms; then we get 



+ p+V(x) = ±u z + p a 



(12) 



Integration of eq. (10 1 over space shows that stationary 



patterns satisfy the condition (see also [14]) 
PiPo - P)&x = 0. 



(13) 



Excluding u from Eqs. (12) and (10 1 we get the equation 



p^ul + 2( Po - p-V{x))\ =2Tp{p G -p) (14) 



which determines the dependence of p on x for a given po- 
tential V(x) provided the solution satisfies the condition 
([9]). This equation represents the so-called hydraulic ap- 
proximation for our system. (It generalizes the well-known 
Thomas- Fermi approximation to non-zero flow velocity.) 
At the tails of the wave pattern (|x| — > 00) we can ne- 



glect the potential V{x) and linearize eq. ( 14 1 with respect 
to small deviations p—po from the steady state. This gives 
the asymptotic behavior 



/ 2rp w 
Po\ oc exp wx 



(15) 



which shows that such tails extending beyond the range 
of the potential can exist only under the condition 



\Po ~ u o\ 
Tp u 



»Z. 



(16) 



This means that the sound velocity c s = ^/po separates 
two different regimes — subsonic (uq < c s ) and supersonic 
(uq > c s ) — with drastically different properties. We shall 
discuss them separately. 

Subsonic flow. — The hydraulic approximation de- 
scribes the profiles of the disturbance well enough practi- 
cally in the entire subsonic regime uq < y/po if the velocity 
is not too close to the sound velocity. In fig. [2] we com- 
pare the direct numerical solution of eq. ([!]) with its hy- 
draulic approximation obtained by the numerical solution 



of eq. ( 14 ) ; quite good agreement is observed (notice that 
in all simulations we use 7 = T = 0.05 that gives sound 
velocity c s = 1). Naturally, the asymptotic behavior |l5| ) 
is also confirmed in the upstream flow x < 0. However, 
there is no linearized solution decaying to zero at x —¥ +00, 
hence the transition to the asymptotic steady flow in the 
downstream region x > has to occur within the range 
of the potential V(x), and this conclusion is confirmed by 
direct numerical solution. Indeed, as one observes in fig. [2] 
the stationary disturbance has a smooth shape with a very 
long monotonically decaying left tail. This tail is attached 
to the density dip located almost completely within the 
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Fig. 2: (Color online) Comparison of density profile at uq — 
0.6, Vm = 0.1, and 7 = Y — 0.05 obtained numerically from 
eq. |TJ (black curve) and in the framework of the hydraulic 
approximation (red curve). The dashed line shows plane wave 
solution in absense of potential. 



potential region. The amplitude of the disturbance (i.e. 
the difference between maximal and minimal density in 
the condensate) in the subsonic regime monotonically in- 
creases with growth of the potential strength. It is conve- 
nient to introduce the "number of particles" 



N = [p 1/2 (x) -a] 2 dx 



(17) 



disturbed by the flow. This variable monotonically in- 
creases with the strength of potential [fig. |3ja)]. Inter- 
estingly, there exists a critical strength of potential V m r 
above which one cannot find stationary profiles. When 
V m — > the tangent to N(V m ) becomes vertical. For 
this strength of potential the amplitude of disturbance 
is maximal. Note that the critical potential strength 
V^ 7 " diverges when the incident velocity of the condensate 
Mo — > and it monotonically decreases when uq increases 
[fig. [3^b)]. We were able to obtain the critical potential 
strength only for < uo < 0.8c s . For higher velocities 
[0.8c s < u < c s ] the condensate starts developing small 
oscillations on its right tail. Solutions characterized by dif- 
ferent number of oscillations on their right tail may form 
continuous families that may be obtained even for poten- 
tials with V m > 1. 



The condition ( 16 ) breaks down when the incident ve- 
in this case the 



locity uo is close to the critical value ^/po 
potential V(x) cannot be neglected in the hydraulic ap- 
proximation ( 14 Even more, as Uq approaches the critical 



value, the upstream density profile steepens and its slope 
can become so large that the dispersive terms in eq. (Ill 



can no longer be neglected. As is known, the dispersive ef- 
fects lead to generation of oscillations in the regions with 
fast change of the variables. Hence, for u > -y/po we 
should expect generation of dispersive shock waves and 
now we proceed to the discussion of the supersonic flow. 



Supersonic flow. 



If uq > Jp~Q and the condition 



(16 1 is fulfilled, then the linearized solution (15) describes 
decaying disturbance in the downstream flow at x — > +00. 
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Fig. 3: (a) Number of particles (17 1 versus potential strength 
for subsonic flow with uo = 0.6. The circle corresponds to the 
parameters used in fig. [5] (b) Critical strength of potential 
versus uq. In all cases 7 = F — 0.05. 



Thus, contrarily to the subsonic case, the right tail for 
x > can extend far beyond the range of the potential. 

As mentioned above, the generation of an oscillatory 
dispersive shock wave is expected in the upstream flow. 
This is confirmed by numerical solution of ([I]) with the 
supersonic boundary conditions; see fig. [4] The dispersive 
shock waves can be represented as a modulated solution 
of the undisturbed (r = 7 = 0, V(x) = 0) equation ([I]) or 
the system Q , Q which can be written in the form (see, 
e.g., M) 



p(x,t) = v% + (v 2 - ^i)sn (V^3 — v\ 9,m) 

= i(A x - A 2 - A 3 + A4) 2 + (A 2 - Ai)(A 4 - A 3 ) 



x sn 2 (v/(A 4 - A 2 )(A 3 - Ai) 9, m), 



(18) 



u{x,t) = U - q/p(x,t), 

g = -§(Ai - A 2 - A 3 + A 4 )(Ai - A 2 + A 3 - A 4 ) 



x (Ai + A 2 - A 3 - A 4 ), 



(19) 



where 



x-Ut, U= i(Ai + A 2 + A 3 + A 4 ) (20) 



v \ ? v 2 , are related to the Riemann invariants Ai , i 
1, 2, 3, 4, by the formulae 



v\ = |(Ai - A 2 - A 3 + A 4 ) 2 , 
»2 = j(Ai -A 2 + A 3 -A 4 ) 2 , 



v-.\ 



— + A 2 — A3 — A 4 



and 



(A 2 -Ai)(A 4 - A 3 ) 



(21) 



(22) 



v z -v x (A 3 - Ai)(A 4 - A 2 ) 

In the modulated wave the Riemann invariants A^ become 
slow functions of x and t and their evolution is governed by 
the Whitham equations. They were derived for the case 
of eq. |T]) with r = in [l7j and their generalization to 
the nonzero T is straightforward. Therefore we will write 
down here the final result without its derivation. 
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Fig. 4: (Color online) The density distribution for shock wave 
at uo = 1.4, V m = 0.8, and 7 = F = 0.05. Red curves show ana- 



lytical approximations of tails from eqs. ( 15 \ and (42 1. Dashed 



line shows plane wave solution in the absence of potential. 

In our stationary case Ai do not depend on t and the 
phase velocity U equals to zero, 

Ai + A 2 + A 3 + A 4 = 0. (23) 

Then the Whitham equations can be written in the form 



hK + h 



where 



dAi _ 2 



(24) 



h = 



Pq-v 



(25) 



H={v- vi){v - v 2 ){v - 1/3) 
and L is the wavelength 

_ 2K(to) 
~ v / (^3-Ai)(A 4 = I 2 7' 

K(m) being the complete elliptic integral of the first kind 



(26) 



(27) 



Due to the special structure of Eqs. (24 1, the symmetric 
functions of Ai , 

s 1 = 2jA i , s 2 = J] XjXj, s 3 = 2J ^i^j^k, 

ijtj i^jjtk (28) 

s 4 = AiA 2 A3A 4 , 



obey very simple equations 




dfi _ ds2 
dx dx 



ds 3 
dx 



2/i ds 4 
L ' da; 



2h 
L ' 



(29) 



The condensate density p oscillates in the interval vi < 
p < ^2, i-e. v\ and v 2 are the envelopes of the oscillations 
in the dispersive shock wave. Let us study their asymp- 
totic behavior at x — > —00. 
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The asymptotic plane wave corresponds to v 2 = v\ (or 
A2 — Ai) so that the parameters of the incident condensate 
are expressed in terms of the Riemann invariants as 



Po 



|(A 4 -A 3 ) 5 



and 



u = i(A 3 + A 4 -2Ai). 



(30) 



(31) 



From (23), (30) and (31) we find the asymptotic values of 



the Riemann invariants 



Xi — X 2 — — 

M — ~z 



2 ' 
Pa, 



. "0 

Y 



(32) 



Po- 



The periodic solution corresponds to the ordering Ai < 
A2 < A3 < A4, so that from A3 > A2 we find that uq > ^fpo, 
i.e. the incident flow must be supersonic for the generation 



of such a stationary shock wave. The wavelength (27) of 



small amplitude oscillations (m = 0) reduces in the limit 
x —00 to 

L= * (33) 

V u o - Po 

The asymptotic values of ^1,^21 v 3 a t x — > —00 are 
V\{— 00) = v 2 {—oo) = po, Us(-oo) = Uq. We introduce 
small deviations from these values, 



P0 + S1/1, v 2 =p + 8v 2 , vi 



+ Sva, (34) 



and linearize eqs. (29) with respect to 8vi with account of 
the identities 



v\ + v 2 + v-3 = -2s 2 , 

V\V% + V1U3 + U2V3 = ~ s \ + 4S4) V\V 2 V3 

As a result we arrive at the equations 
d(Svi + 8v 2 ) d<5^3 ZTpquq 



(35) 



dx 

Hence we get 



dx 



Po 



{8v 1 +8v 2 ). (36) 



oz/i + 0^2, 0^3 oc exp — ^ x 



(37) 



These equations suggest that the sum of deviations 8i>i + 
8v 2 decays at x — ¥ — 00 faster than the deviations 8vi and 
8v 2 separately. If we introduce small deviations of the 



Riemann invariants from their asymptotic values ( 32 ) 
Ai 



' 2 



+ 8\i, \ 2 



uo 
2 



8\o 



A3 = -y - Vpo + 8X3, A 4 = -y + y/po~ + 8X 4 , 



(38) 



then we find that the asymptotic behavior (37) corre- 
sponds to 1 5 A3 1 , \5\4\ <C \8Xi\, \8X 2 \ so that 



i(2Vno + ^Ai-5A 2 ) 5 



»i = \{2^ ~ 8X 1 + 8X 2 f 



(39) 



v z = i(2 Uo -JA 1 -5A 2 ) 2 




Fig. 5: (a) Number of particles versus potential strength for 
supersonic flow with ito = 1.4. The circle corresponds to solu- 
tion in fig. [4] (b) Critical strength of potential versus no. In 
all cases 7 = Y = 0.05. 



and in leading order approximation 8v\ + 8v 2 = which 



means that (37) corresponds to higher order corrections. 



Thus, in the main approximation we assume that 2/3 = 



Uq and from (25 1 we obtain the following equations for <5Ai 
and 8X 2 : 



d<5Ai 



Tp u Q 



dx 2{ul - po) 
d8X 2 _ Tp u 
dx 2(ul - p a 

From these equations we get 
8\\ 2 = ±C exp 



{8X 1 - 5X 2 ), 
(8X 1 -8X 2 ). 



(40) 



TpoUo_ 
«o - Po ' 



(41) 



(C is the integration constant), and hence the envelopes of 
oscillations in the dispersive shock decay at leading order 
approximation as 



8v\, 8v 2 oc exp 



(42) 



This is slower than the dependence in ( 37 1 as it should be 



These analytical predictions were confirmed by the di- 
rect solution of eq. ([I]). As one can see from fig.|4j a typical 
shock wave solution may extend far beyond the repulsive 
potential, both in the positive and negative x directions. 
Note the excellent agreement at x — > ±00 between the 
numerical density and its analytical fit that uses the ex- 
pressions (15) and (42). The wavelength of small ampli- 



tude oscillations estimated numerically for the parameters 
of the flow shown in fi g. |4| is equal to L num = 3.16 and 
the analytical formula (33) gives L t heor — 3.21, with an 
accuracy better than 2%. It should be stressed that the 
solution in fig. [4] is fully stationary (that is, it does not 
change in time), in contrast to previously reported one- 
dimensional shock waves in conservative systems. The os- 
cillating left tail of the dispersive shock wave in the super- 
sonic regime becomes more pronounced with increase of 
potential strength V m and incident velocity uq. Just as for 
the subsonic flow the increase of V m is accompanied by an 
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increase of the amplitude of the shock wave. However, in 
the supersonic regime the amplitude of shock wave may be 
comparable with the amplitude of the unperturbed plane 
wave and the density may decrease almost to zero in the 
downstream region, especially for strong potentials. The 
number of particles in the shock wave increases with V m 
[fig. (5ja)] and the character of this dependence also points 
to the existence of a critical defect strength beyond which 
shock waves do not exist. 

A linear stability analysis performed for stationary solu- 
tions shows that wave patterns are stable for any strength 
of potential V m up to the critical one as long as the incident 
velocity is small enough, u < 0.9c s . When u > 0.9c s , 
stationary solutions become unstable if the strength of po- 
tential exceeds certain limiting value V£[ (the instability of 
the dispersive shock waves at V m > is accompanied by 
the development and emission of small-scale disturbances 
on the right tail of the shock wave). Still, even shock waves 
with very long oscillating tails, like the one shown in fig|4j 
may be stable. The critical value of potential increases 
monotonically with uq as it is shown in fig.[5jb). 

Conclusion. — It is instructive to compare our results 
with those obtained in the case of the flow of a conserva- 
tive fluid past an obstacle. As was shown in 



12 , in the 



conservative case one can distinguish three characteristic 
ranges of flow velocity — subcritical (uq < transcrit- 
ical (u_ < uq < and supercritical (uq > u+) where 
the critical values of the velocity are located at opposite 
sides of the value of the sound velocity (m_ < c s < u + ). 
If the flow velocity is subcritical or supercritical, then the 
disturbance is stationary and has the dimension of the ob- 
stacle's size with definite sign of the difference p—po; it has 
a form of a dip in subcritical region and of a hump in su- 
percritical region. In the present case of non- conservative 
dynamics the disturbance must have "dips" and "humps" 
to be compensated in the integral ( [13] ) . Examples of such 
disturbances shown in figs. [2] and [4] demonstrate this prop- 
erty. 

In the conservative situation the transcritical regime 
corresponds to non-stationary generation of upstream 
and/or downstream dispersive shock waves. In the non- 
conservative situation this regime actually disappears and 
instead a supersonic flow generates a stationary dispersive 
shock wave upstream the obstacle. There are no down- 
stream dispersive shock waves; instead they are replaced 
by smooth profiles decaying asymptotically to the station- 
ary plane wave. 

The above properties of the wave patterns generated by 
the flow of polariton condensate past an obstacle are de- 
rived for the model ([I]). We believe, however, that these 
properties remain qualitatively the same for other models 
having stationary plane wave states formed due to balance 
of pumping and dissipation effects. In particular, this the- 
ory can find applications to a so-called "superfluid motion 
of light" fiMM 
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